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Many-body localization and thermalization in disordered Hubbard chains 
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We study the many-body localization transition in one-dimensional Hubbard chains using exact 
diagonalization and quantum chaos indicators. We also study dynamics in the delocalized (ergodic) 
and localized phases and discuss thermalization and eigenstate thermalization, or the lack thereof, 
in such systems. Consistently within the indicators and observables studied, we find that ergodicity 
is very robust against disorder, namely, even in the presence of weak Hubbard interactions the 
disorder strength needed for the system to localize is large. We show that this robustness might be 
hidden by finite size effects in experiments with ultracold fermions. 
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Introduction. Over the years, substantial attention 
has been devoted to understanding the dynamical prop¬ 
erties of disordered systems. Interest on this topic goes 
back to a seminal paper by Anderson in 1958, who 
showed that sufficiently strong quenched disorder can 
produce localization of noninteracting particles, preclud¬ 
ing transport in the thermodynamic limit [1]. Destruc¬ 
tive interference is at the heart of this phenomenon. It 
is more prominent in lower dimensions, and, as a result, 
any nonzero disorder strength leads to localization in one 
and two dimensions [2]. A fundamental aspect of Ander¬ 
son localization is that it occurs not only in the ground 
state but also in (highly) excited states. 

Because of the possibility of localization occurring in 
interacting systems, a phenomenon termed many-body 
localization (MBL), disordered systems in the presence 
of interactions have received a lot of attention in recent 
years. Early perturbative arguments [3-6] and numerical 
simulations in the presence of strong interactions [7-11] 
have triggered much research on this topic [12, 13]. The 
MBL transition has also started to be explored in exper¬ 
iments with ultracold atoms [14-16] and ions [17]. 

The contrast between the properties of many-body 
eigenstates of interacting systems in the presence and ab¬ 
sence of MBL makes apparent how remarkable MBL is. 
In generic isolated systems, interactions make it possible 
for the system to act as its own “effective bath”. If taken 
out of equilibrium, such systems evolve in time in such a 
way that observables equilibrate and can be described by 
traditional ensembles of statistical mechanics (i.e., they 
thermalize). This is just one of the manifestations of a 
phenomenon known as eigenstate thermalization [18-20], 
which, in short, means that the expectation value of an 
observable in an eigenstate of a many-body interacting 
system is the same as that in thermal equilibrium (with 
the same mean energy as the eigenstate energy). Eigen¬ 
state thermalization has been shown to occur in several 
many-body quantum systems [20-29]. It is known not to 
occur only in integrable and MBL systems, i.e., the latter 
two classes of systems generally do not exhibit thermal¬ 
ization even if they are thermodynamically large [30, 31]. 

As a matter of fact, it was the latter property of MBL 


systems that was used in the experiments of Ref. [15] 
to distinguish between the delocalized (ergodic) regime 
and the MBL one, for spinful fermions in the presence 
of a quasi-periodic potential. Motivated by those experi¬ 
ments, in this work we study the MBL transition in Hub¬ 
bard chains with disorder. We contrast the predictions of 
quantum chaos indicators for the transition to those from 
thermalization and eigenstate thermalization. We argue 
that ergodicity is remarkably robust in these itinerant 
systems, and show that finite size effects in thermaliza¬ 
tion indicators might hide this fact in experiments. 

Model and the MBL transition. To investigate the 
MBL transition, we use full exact diagonalization and 
study the Hamiltonian: H = Ho + H s b + H\y , in which 
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is an extended Hubbard model (written in standard no¬ 
tation) in a linear chain of size L (with open bound¬ 
ary conditions), with nearest neighbor hoppings (am¬ 
plitude f), onsite interaction (strength U), and next- 
nearest neighbor hoppings (amplitude t'). We have taken 
t' / 0 so that the model is nonintegrable (quantum 
chaotic) in the absence of disorder. Additional sym¬ 
metries, parity, and SU(2), are removed by adding a 
very weak magnetic field (hb) and chemical potential 
(fib), respectively, at the opposite edges of the chain: 
Hsb = M”i,t - ”-1,4.) + Ub(fiL ,t + h L , 4 .) (see Ref. [32] 
for details). We focus on a uniformly distributed dis¬ 
order described by Hw = X)i C r £ *”*c r ’ where the local 
potential gy g [— W/2, W/2]. To contrast the effect of 
the disorder with the effect of the quasi-periodic poten¬ 
tial studied in Refs. [15], we also report results for the 
phase diagram when £» = y cos (2n[3i + </>), where A is 
the potential strength, /? = (y/5 + l)/2 is the golden 
ratio, and <f> is an arbitrary phase (as in the Aubry- 
Andre model [33]). Throughout this Rapid Connnuni- 
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cation, t = 1 sets the energy scale and t' = 0.5. We only 
change U and the disorder strength. The systems stud¬ 
ied are at quarter filling, namely, + Aq = L/2, with 
N t = (J2f hi t ) and N± = We consider two 

lattice sizes, L = 10 and 12, where N-f = = L/4 for 

L = 12, and N t = N± ± 1 for L = 10 [32], 

A common quantum chaos indicator used to locate the 
many-body localization transition in disordered systems 
is the average ratio between the smallest and the largest 
adjacent energy gaps, r n = min [5®, max[<5f, 

with 5% = E n — E r j_i, and {E n } is the ordered list of en¬ 
ergy levels [7]. Here, in order to reduce finite size effects, 
we compute the average ratio f over the central half of the 
spectrum. In the ergodic phase, when the level spacing 
exhibits a Wigner-Dyson distribution, the average ratio 
is r WD « 0.536, while in the MBL phase, when the level 
spacing exhibits a Poisson distribution, the average ratio 
is r p = 2In2 — 1 « 0.386 [34], 

Figure 1 shows the disorder average of f, (f)d; s , as a 
function of the disorder strength for different values of 
the on-site repulsion and two system sizes. The value of 
the disorder strength at which the curves cross or merge, 
W c , can be taken as an estimate of the critical disorder 
strength for the ergodic to many-body localization phase 
transition. Such a crossing or merging point is known to 
move towards stronger disorder with increasing system 
size (see, e.g., Refs. [7, 9, 10]); as such, the values re¬ 
ported here should be thought of as lower bounds for the 
critical disorder. As expected, since interactions promote 
delocalization, W c first increases with U [Figs. l(a)-l(c)]. 
It is remarkable that, even for fairly small values of U 
[U = 0.2 in Fig. 1(a)], the delocalized regime is robust 
up to values of W c ~ 8, i.e., almost twice the width B of 
the single particle spectrum, e*, = — 2f cos(fc) — 2f' cos(2fc) 
(B = 4.5 for our parameters). As the onsite interaction 
strength becomes of the order of B , W c stops increas¬ 
ing and, as U increases further, W c starts to decrease 
[Figs. 1(c) and 1(d)]. This is expected as, in the limit 
U —¥ oo, each sector in the Hubbard model with a par¬ 
ticular ordering of the spins (and no double occupancy) 
maps onto a noninteracting spinless fermion Hamiltonian 
with Nf + Ni fermions, and the latter localizes for any 
nonzero disorder strength. Figure 1(e) depicts the esti¬ 
mated phase diagram in the presence of disorder for up 
to U = 20. In contrast, as also shown in Fig. 1(e), MBL 
in the presence of a quasi-periodic potential (see Ref. [32] 
for further details) occurs for A < W. MBL is also easier 
to achieve in interacting spinless fermion systems [35]. 

Dynamics and thermalization. As mentioned before, 
one of the defining properties of the MBL phase is its 
lack of thermalization. In what follows, motivated by 
the experimental results reported in Ref. [15], we study 
dynamics in the delocalized and MBL regimes. Our ini¬ 
tial state is also experimentally motivated. We consider 
| ifii) = | f 0 4- 0 t 0 4- ■••); which is a state that has 
no double occupancy and can be prepared using optical 



FIG. 1. (Color online) (a)-(d) Averaged ratio of adjacent 
energy gaps as a function of the disorder strength for four 
values of U and two lattice sizes. The average f was com¬ 
puted over the central half of the spectrum. The disordered 
averaged results (f)dis for L = 10 were obtained averaging 
over 1200 disorder realizations, and the ones for L = 12 over 
20-200 disorder realizations (error-bars report the standard 
deviation). In (a), we show results for Nf = -W ± 1 when 
L = 10. They make apparent that both sectors behave qual¬ 
itatively (and quantitatively) similarly even for the largest 
values of hb = fib used. The crossing (or merging) point be¬ 
tween curves for different lattice sizes provides an estimate of 
the critical disorder, W c , for the ergodic to MBL transition, 
(e) Estimated W c and A c as a function of U (error bars report 
an interval of confidence based on the closeness of the results 
for L = 10 and L = 12 about W c and A c ). 


superlattices. | ipi) is a quarter-filling version of the state 
prepared in Ref. [15]. The dynamics is then studied un¬ 
der H = H 0 + I7 s b + Hw- Our goal is to understand 
how the results of the dynamics relate to those obtained 
for f. Some of the specific questions we address are the 
following: Is the MBL transition manifest in the dynam¬ 
ics of experimentally relevant observables? At what time 
do those observables reach (if they do) stationary val¬ 
ues? We are also interested in understanding the role of 
finite size effects. They have been found to be stronger 
in indicators related to Hamiltonian eigenstates than in 
those related to the spectrum [36]. To address these ques¬ 
tions, we focus on one particular value of the interaction 
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FIG. 2. (Color online) Disorder averaged results for the time 
evolution of (a) the even-odd site occupation imbalance, and 
(b) the antiferromagnetic structure factor. The shaded area 
around the curves depicts the standard deviation of the mean, 
after an average over ten disorder realizations. The horizon¬ 
tal dashed lines depict the disorder averaged values of the 
diagonal ensemble predictions (see text). 


strength, U = 4. 

We report results for three observables (see Ref. [32] for 
another one). Two observables, the imbalance I = ((h e )— 
(n°»/((n e ) + <n°)), where h e = Ei= e ve„(odd),<7 ( I 

was measured in Ref. [15]), and the kinetic energy K = 

+ H.C.)) - a Ei, CT ((Cci+2, t 7 + H.c.)), 
are directly related to the charge degrees of freedom. 
The third one, the antiferromagnetic structure factor 
S = l/fE,j e ' l( '" : ' ) ((«it - ~ Rr)) is related 

to the spin degrees of freedom (from now on we refer to 
it as the structure factor). The relaxation times of the 
charge and spin degrees of freedom are expected to be 
different for very strong interactions [37]. 

Figure 2 displays the disorder averaged time evo¬ 
lution of the imbalance [Fig. 2(a)] and of the struc¬ 
ture factor [Fig. 2(b)]. We also display, as horizontal 
dashed lines, the disorder average of the diagonal en¬ 
semble results. Given an observable O, the diagonal en¬ 
semble result (which, in the absence of degeneracies, is 
the same as the infinite-time average of the observable 


FIG. 3. (Color online) (a) Disorder averaged diagonal en¬ 
semble results for the imbalance (main panel) and the kinetic 
energy per site (inset) vs the amplitude of the disorder W. 
(b) Normalized disorder average difference between the diago¬ 
nal and microcanonical ensemble predictions for the structure 
factor (main panel) and the kinetic energy (inset) vs the am¬ 
plitude of the disorder W. In all cases U = 4, and the width 
of the microcanonical energy window is A E = 0.1. The ver¬ 
tical dashed line marks W c , and the shaded region around it 
signals the interval of confidence reported in Fig. 1(e). 


[20]) can be obtained as O de = E Q |C , a | 2 C , Q , a , where 
O aa = (a|0|a), | a) are the eigenstates of the Hamil¬ 
tonian (H |a) = E a |a), E a are the eigenenergies), and 
C a = (a\xfi). We say that O equilibrates if it relaxes to 
Ode and remains close to it at later times. 

Figures 2(a) and 2(b) show that I and S equilibrate in 
both the delocalized and localized regimes. Sufficiently 
far away from W c (W = 0, 2, and 30 in the figure), one 
can see that both observables have essentially reached 
the diagonal ensemble result (or are very close to it) 
for t ~ 10 (h/t). As the system approaches W c (~ 20 
for U = 4), we find that equilibration times become 
much longer. For example, for W = 16 in Fig. 1(a), one 
can see that I becomes nearly time independent only at 
r ~ 10 3 (h/t). Very long equilibration times at a delocal¬ 
ization to localization transition have also been observed, 
for much larger system sizes, in the integrable hard-core 
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boson version of the Aubry-Andre model [38]. Those 
times represent a challenge for experiments. 

Next, we check how the diagonal ensemble results for 
the observables compare to the microcanonical predic¬ 
tions. Whenever they agree, and equilibration occurs (as 
we have checked), we say that the system thermalizes. 
We first consider I. Since there is no distinction between 
even and odd sites in the Hamiltonian, the disorder av¬ 
erage of / is expected to be zero in the microcanonical 
ensemble ((/me) dis = 0)- Hence, as argued in Ref. [15], 
the disorder average of the diagonal ensemble result for 

1 ((Joe) dis) can be taken to be the order parameter for 
the MBL phase (it can only differ from zero if the system 
does not thermalize). Figure 3(a) shows (iDE)dis vs W 
for two system sizes. For the (small) system sizes that we 
can study, (/oE)dis can be seen to smoothly increase from 
zero with increasing W. However, comparing the results 
for the two system sizes, one can see that in the delo¬ 
calized side (and close to W c in the MBL side) (/DE)dis 
decreases with increasing system size. This is consistent 
with the expectation that, in the thermodynamic limit, 
it will vanish in the delocalized side. 

Another order parameter that could be used to locate 
the MBL transition in experiments is the kinetic energy. 
As discussed in Ref. [31], the dynamics of one-particle 
correlations in the MBL phase is quantitatively similar to 
that in the atomic limit (even if the system is not close to 
that limit). This means that, in the Heisenberg represen¬ 
tation, c| i(7 (r)c ji a (r) « exp[i(£i - e J )r/fi]cJ )(T (0)c jiCT (0). 
Given our initial state, that implies that (A’DE)dis ~ 0. 
In the inset in Fig. 3(a), one can see that, indeed, 
(ATDE)dis ~ 0 for W > W c . In Fig. 3(b), the inset shows 
(|A"de - A*me|) dis/(|A" me|) dis and the main panel shows 
(|S'DE-*S , ME|)dis/(|-S , ME|)dis- Both normalized differences 
can be seen to decrease with increasing system size in the 
delocalized phase and not in the MBL. This is consistent 
with the expectation that thermalization occurs only in 
the former. 

The fact that the system thermalizes (fails to thermal¬ 
ize) in the delocalized (MBL) regime can be understood 
to be the result of eigenstate thermalization occurring 
(not occurring) in that regime [18-20]. In Fig. 4, we show 
the eigenstate expectation values of the three observables 
of interest for a single disorder realization for different 
values of W. Deep in the delocalized phase (W = 0 and 

2 in the figure), the support for those expectation val¬ 
ues at a given energy can be seen to be very small (it 
decreases with system size, not shown), i.e., eigenstate 
thermalization occurs. The support of the eigenstate ex¬ 
pectation values exhibits a different behavior within the 
MBL phase, or close to it (W = 16 and 30 in the figure), 
i.e., eigenstate thermalization does not occur, or at least, 
it is not apparent for the system sizes studied. In Fig. 4, 
vertical dashed lines depict the mean energy, and shaded 
areas around them depict the width of the energy dis¬ 
tribution (for W = 0 and 2), in the quenches involving 



E a /(U+W) 


FIG. 4. (Color online) Eigenstate expectation values of the 
imbalance (a), the kinetic energy per site (b), and the struc¬ 
ture factor (c) for a single disorder realization in systems with 
four different disorder strengths. The vertical dashed lines 
show the averaged mean energy, and the shaded are for W = 0 
and 2 and depict the averaged energy width, for the quenches 
involving this disorder realization. The results reported here 
were obtained in systems with U = 4 and L = 12. 


that disorder realization. They show which part of the 
spectrum is relevant to the dynamics studied. 

Summary and discussion. We have studied the er- 
godic to MBL transition in Hubbard chains. Our main 
result from the analysis of quantum chaos indicators is 
that ergodicity is very robust against disorder. Even for 
on-site interactions as weak as U = 0.2, we find that the 
disorder strength required to localize the system is or the 
order of twice the single-particle bandwidth. We studied 
the dynamics of those systems starting from a state of 
the form | ipi) = | fO^OtOJ....). We find that various 
experimentally relevant observables equilibrate in time 
scales ~ 1 — 10 (h/t), whenever the system is not close to 
the MBL transition. Close to the MBL transition, equi¬ 
libration times become orders of magnitude longer and 
might be difficult to reach experimentally. We have also 
studied the differences between observables after equi¬ 
libration and the predictions of the microcanonical en¬ 
semble finding that, for the small lattice sizes that we 
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are able to study (~ 12 sites), they increase smoothly 
as one increases disorder and can be large even far from 
the MBL transition. This is reminiscent of the behavior 
observed as one approaches an integrable point in finite 
systems [21, 22]. Hence, the analysis of a few small sys¬ 
tem sizes does not allow one to identify the critical disor¬ 
der strength at which MBL occurs. Large system sizes, 
or a careful finite size scaling analysis, are needed. While 
that might be possible in experiments, it remains a chal¬ 
lenge for numerical simulations. Numerical linked cluster 
expansions [30, 31], which exhibit an exponentially fast 
convergence with increasing the size of the systems that 
need to be diagonalized [39], offer a promising way to 
address this challenge [31]. 
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Supplementary Materials: 

Many-body localization and thermalization in disordered Hubbard chains 


LATTICE SIZES AND THEIR HILBERT SPACE 

In Table SI, we show the Hilbert spaces for the differ¬ 
ent lattice sizes that we can study at quarter filling and 
at half-filling. At quarter-filling, the Hilbert spaces for 
L = 6 and L = 8 are very small and, as shown in the 
next section, exhibit erratic behavior. Therefore, in the 
main text we have focused on the lattices with L = 10 
and 12. 

TABLE SI. Lattice configurations and size of the Hilbert 
space ( D ) used in this work. 


L 

71 / 

Nf-Ni 

D 

6 

3 

±1 

90 

8 

4 

0 

784 

10 

5 

±1 

5400 

12 

6 

0 

48400 

8 

8 

0 

4900 

10 

10 

0 

63504 


SYMMETRY-BREAKING FIELD 

While the extended Hubbard model in Eq. (1) is non- 
integrable for nonzero t 1 , in our lattice geometry (which 
has open boundary conditions) there are still two sym¬ 
metries, parity and 517(2) when Nf = N±, that need to 
be removed for the level spacing in the system to exhibit 
a Wigner-Dyson distribution. This is achieved using the 
symmetry breaking terms described under H s b. 


Figure SI shows the dependence of the average ratio 
of adjacent energy gaps for two values of U and different 
lattice sizes (in the absence of disorder). For simplic¬ 
ity, we have selected hb = /ib, and, as in the main text, 
t' = 0.5. Also, to avoid larger finite size effects related 
to the low and high energy edges of the spectrum, the 
averages are calculated over the central half of the spec¬ 
trum. Due to the smallness of their Hilbert spaces, the 
systems with L = 6 and 8 exhibit very large fluctua¬ 
tions in Fig. SI (this is the reason they were not shown 
in the main text). For very small and very large values 
of hb = /it,, one can see that the average ratio of ad¬ 
jacent energy gaps is consistent with that of a Poisson 



FIG. SI. (Color online) Average ratio of adjacent energy gaps 
as a function of the amplitude of the symmetry breaking terms 
hb = h'b for a chain with U = 1 (a) and [7 = 4 (b), in the 
absence of disorder. In both cases t' = 0.5 and the average is 
computed over the central half of the spectrum. 
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distribution. For small values of hb = fit,, the reason are 
the symmetries mentioned above while, for large values 
of hi, = fib, the reason is that the sites with the fields 
decouple from the rest of the chain. It is apparent in 
Fig. SI that, as the system size increases, the strength 
of the symmetry breaking terms required to produce a 
distribution of level spacings that is consistent with the 
Wigner-Dyson distribution decreases. Hence, for ther¬ 
modynamically large systems, arbitrary weak symmetry 
breaking fields would be required. 

For the results reported in the main text, we use values 
of hb,fib in the range [0.1, 0.5]. The largest values of the 
breaking field were necessary in the limit of weak on-site 
interactions, as the system approached the noninteract¬ 
ing (integrable) limit. 


QUASI-PERIODIC POTENTIAL 

We have also studied the MBL transition in the 
presence of a quasi-periodic potential with £j = 
jy cos (2n/3i + (f>), where A is the potential strength, /3 = 
(-\/5 + l)/2 is the golden ratio, and (j) is an arbitrary phase 
that we use to generate different realizations of {e^} (we 



FIG. S2. (Color online) Average ratio of adjacent energy gaps 
as a function of the strength of the quasi-periodic potential 
A, for U = 4, in systems at quarter filling (a) and half filling 
(b). In both cases, t' = 0.5, f is computed over the central 
half of the spectrum, and (•)</, is an average computed over 
random phases cj>. 


take (j) £ [—7r, 7r] with a uniform distribution). For this 
quasi-periodic potential, unlike for the Anderson case, 
localization in the noninteracting limit occurs only for 
A > 4 (in units of t = 1). As for the disordered case, 
the phase diagram for this potential reported in Fig. 1(e) 
was obtained via an analysis of the average ratio of ad¬ 
jacent energy gaps (see Fig. S2). This phase diagram 
was explored experimentally in a recent optical lattice 
experiment [15]. 

Figure S2 reports the equivalent of Fig. l(a)-l(d) in 
the main text but for the quasi-periodic potential when 
the system is at quarter filling [(a)—(d)] and at half-filling 
[(e)-(h)]. In the latter case, the value of A required to 
drive the MBL transition is slightly larger in comparison 
to the former case. This is understandable as U plays 
an increasingly important role with increasing the filling. 
We expect a similar behavior to occur in the presence 
of disorder. Namely, W c at half-filling should be larger 
than at quarter filling, which was the case discussed in 
the main text. 


OTHER QUANTITIES SUPPORTING MBL 


Other quantum chaos indicators that can be used to 
locate the ergodic to MBL transition are the Shannon 
entropy, also known as the information entropy, 

S a = ~ ^El4l 2l n|4l 2 j > (2) 


and the inverse participation ratio (IPR), 


ipr q 


1 



(3) 


where D is the dimension of the Hilbert space, and the 
coefficients c J a correspond to the _y-t.li component of the 
energy eigenvector |a) in some basis. Here, we use the 
computational (Fock) basis. As one increases the system 
size in the ergodic (chaotic) regime, one expects S a to 
approach the Gaussian orthogonal ensemble (GOE) pre¬ 
diction S2° E = ln(0.48D), provided |a) is away from the 
edges of the spectrum [36]. For IPR a one expects that it 
should approach IPR GOE = (D + 2)/3. Figure S3 shows 
that, when W = 0 and 2 for L = 12, both S a and IPR a 
are indeed close to the GOE predictions away from the 
edges of the spectrum. For W = 16 and 30, on the other 
hand, localization in the Fock basis is made apparent by 
the small values of S a and IPR a over the entire spectrum. 

In Fig. S4, we show the energy eigenstate expectation 
values of another observable that is of interest to experi¬ 
ments with ultracold fermions in optical lattices, namely, 
the double occupancy: n-|q = 1/L The re¬ 

sults for this observable are qualitatively similar to those 
discussed in the main text (see Fig. 4 there) for the im¬ 
balance, the kinetic energy, and the structure factor. 
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FIG. S3. (Color online) (a) Shannon entropy of the eigen¬ 
states of the Hamiltonian in the computational (Fock) basis 
vs the eigenstate energies, (b) Same as (a) for the inverse 
participation ratio. The parameters are: L = 12, U = 4, and 
t' = 0.5. Results are presented for a single realization of dis¬ 
order and four values of W. The horizontal lines depict the 
GOE predictions. 
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FIG. S4. (Color online) Same as Fig. 4 in the main text but 
for the double occupancy. 















